home *** CD-ROM | disk | FTP | other *** search
/ Aminet 5 / Aminet 5 - March 1995.iso / Aminet / dev / misc / LEDA_gene.lha / LEDA-3.1c-generic / prog / plane / seg_int_rat.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-05  |  13.8 KB  |  616 lines
  1. #include <LEDA/sweep_segments.h>
  2. #include <LEDA/sortseq.h>
  3. #include <LEDA/prio.h>
  4. #include <math.h>
  5. #include <LEDA/rational.h>
  6.  
  7.  
  8. /*
  9. #define EPS  0.00001
  10. #define EPS2 0.0000000001
  11. */
  12.  
  13. #define double rational
  14.  
  15. const double EPS  = 0;
  16. const double EPS2 = 0;
  17.  
  18. class SW_POINT;
  19. class SW_SEGMENT;
  20. typedef SW_POINT* SW_point;
  21. typedef SW_SEGMENT* SW_segment;
  22.  
  23. enum SW_point_type {Cross=0,Rightend=1,Leftend=2}; 
  24.  
  25. class SW_POINT
  26. {
  27.   friend class SW_SEGMENT;
  28.  
  29.   SW_segment seg;
  30.   int     kind;
  31.   double  x;
  32.   double  y;
  33.  
  34.   public:
  35.  
  36.   SW_POINT(double a,double b)  
  37.   { 
  38.     x=a; y=b; seg=0; kind=Cross;
  39.    }
  40.  
  41.  
  42.  
  43.   SW_POINT(point p)         
  44.   { 
  45.     x=p.xcoord();y=p.ycoord();seg=0;kind=Cross;
  46.    }
  47.  
  48.  
  49.   LEDA_MEMORY(SW_POINT);
  50.  
  51.   friend double    get_x(SW_point);
  52.   friend double    get_y(SW_point);
  53.   friend int       get_kind(SW_point);
  54.   friend SW_segment get_seg(SW_point);
  55.  
  56.   friend bool intersection(SW_segment, SW_segment, SW_point&);
  57. };
  58.  
  59. inline double    get_x(SW_point p)    { return p->x; }
  60. inline double    get_y(SW_point p)    { return p->y; }
  61. inline int       get_kind(SW_point p) { return p->kind; }
  62. inline SW_segment get_seg(SW_point p)  { return p->seg; }   
  63.  
  64. static int compare(SW_point& p1,SW_point& p2)
  65. { if (p1==p2) return 0;
  66.  
  67.   int c = compare(get_x(p1), get_x(p2));
  68.  
  69.   if (c) return c;
  70.  
  71.   int diffk = get_kind(p1)-get_kind(p2);
  72.  
  73.   if (diffk != 0) return diffk;
  74.  
  75.   return compare(get_y(p1), get_y(p2));
  76.  
  77. }
  78.  
  79.  
  80.  
  81. /*
  82.  
  83. static int compare(SW_point& p1,SW_point& p2)
  84. { if (p1==p2) return 0;
  85.  
  86.   double diffx = get_x(p1) - get_x(p2);
  87.   if (diffx >  EPS2) return  1;
  88.   if (diffx < -EPS2) return -1;
  89.  
  90.   int    diffk = get_kind(p1)-get_kind(p2);
  91.   if (diffk != 0) return diffk;
  92.  
  93.   double diffy = get_y(p1) - get_y(p2);
  94.   if (diffy >  EPS2) return  1;
  95.   if (diffy < -EPS2) return -1;
  96.  
  97.   return 0;
  98.  
  99. }
  100. */
  101.  
  102.  
  103. class SW_SEGMENT
  104. {
  105.   SW_point startpoint;
  106.   SW_point endpoint;
  107.   double  slope;
  108.   double  yshift;
  109.   node  left_node;
  110.   int   orient;
  111.   int   color;
  112.   int   name;
  113.  
  114.  
  115.   public:
  116.  
  117.   SW_SEGMENT(SW_point, SW_point,int,int);     
  118.  
  119.  ~SW_SEGMENT() { delete startpoint; delete endpoint; }     
  120.  
  121.   LEDA_MEMORY(SW_SEGMENT);
  122.  
  123.   friend SW_point get_startpoint(SW_segment);
  124.   friend SW_point get_endpoint(SW_segment);
  125.   friend double get_slope(SW_segment);
  126.   friend double get_yshift(SW_segment);
  127.   friend int  get_orient(SW_segment);
  128.   friend int  get_color(SW_segment);
  129.   friend int  get_name(SW_segment);
  130.   friend node get_left_node(SW_segment);
  131.   friend void set_left_node(SW_segment, node);
  132.  
  133.   friend bool intersection(SW_segment, SW_segment, SW_point&);
  134. };
  135.  
  136.  
  137. inline SW_point get_startpoint(SW_segment seg)     { return seg->startpoint; }
  138. inline SW_point get_endpoint(SW_segment seg)       { return seg->endpoint; }
  139. inline double get_slope(SW_segment seg)           { return seg->slope; }
  140. inline double get_yshift(SW_segment seg)          { return seg->yshift; }
  141. inline int  get_orient(SW_segment seg)            { return seg->orient; }
  142. inline int  get_color(SW_segment seg)             { return seg->color; }
  143. inline int  get_name(SW_segment seg)              { return seg->name; }
  144. inline node get_left_node(SW_segment seg)         { return seg->left_node; }
  145. inline void set_left_node(SW_segment seg, node v) { seg->left_node = v; }
  146.  
  147.  
  148.  
  149. SW_SEGMENT::SW_SEGMENT(SW_point p1,SW_point p2,int c, int n)    
  150.   {
  151.     left_node  = nil;
  152.     color      = c;
  153.     name       = n;
  154.  
  155.     if (compare(p1,p2) < 0)
  156.      { startpoint = p1; 
  157.        endpoint = p2; 
  158.        orient = 0;
  159.       }
  160.     else
  161.      { startpoint = p2; 
  162.        endpoint = p1; 
  163.        orient = 1;
  164.       }
  165.  
  166.     startpoint->kind = Leftend; 
  167.     endpoint->kind = Rightend; 
  168.     startpoint->seg = this; 
  169.     endpoint->seg = this;
  170.  
  171.     if (endpoint->x != startpoint->x)
  172.     {
  173.       slope = (endpoint->y - startpoint->y)/(endpoint->x - startpoint->x);
  174.       yshift = startpoint->y - slope * startpoint->x;
  175.  
  176.       startpoint->x -= EPS;
  177.       startpoint->y -= EPS * slope;
  178.       endpoint->x += EPS;
  179.       endpoint->y += EPS * slope;
  180.     }
  181.     else //vertical segment
  182.     { startpoint->y -= EPS;
  183.       endpoint->y   += EPS;
  184.      }
  185.   }
  186.  
  187.  
  188. static double x_sweep;
  189. static double y_sweep;
  190.  
  191.  
  192. static int compare(SW_segment& s1,SW_segment& s2)
  193. {
  194.   double y1 = get_slope(s1)*x_sweep+get_yshift(s1);
  195.   double y2 = get_slope(s2)*x_sweep+get_yshift(s2);
  196.  
  197.   double diff = y1-y2;
  198.  
  199.   if (diff >  EPS2) return  1;
  200.   if (diff < -EPS2) return -1;
  201.  
  202.   if (get_slope(s1) == get_slope(s2)) 
  203.         return compare(get_x(get_startpoint(s1)), get_x(get_startpoint(s2)));
  204.  
  205.   if (y1 <= y_sweep+EPS2)
  206.         return compare(get_slope(s1),get_slope(s2));
  207.   else
  208.         return compare(get_slope(s2),get_slope(s1));
  209.  
  210. }
  211.  
  212.  
  213. void Print(SW_segment& x) 
  214. { SW_point s = get_startpoint(x);
  215.   SW_point e = get_endpoint(x);
  216.   cout << 
  217.     string("[(%f,%f)----(%f,%f)]",float(get_x(s)),float(get_y(s)),
  218.                                   float(get_x(e)),float(get_y(e)));
  219. }
  220.  
  221.  
  222. static priority_queue<seq_item,SW_point> X_structure;
  223.  
  224. static sortseq<SW_segment,pq_item> Y_structure;
  225.  
  226.  
  227. bool intersection(SW_segment seg1,SW_segment seg2, SW_point& inter)
  228. {
  229.   if (seg1->slope == seg2->slope)
  230.     return false;
  231.   else
  232.   {
  233.     //x-coordinate  of the intersection
  234.     double Cross_x = (seg2->yshift - seg1->yshift) / (seg1->slope - seg2->slope);
  235.  
  236.     if (Cross_x <= x_sweep) return false;
  237.  
  238.     double s1 = seg1->startpoint->x;
  239.     double s2 = seg2->startpoint->x;
  240.     double e1 = seg1->endpoint->x;
  241.     double e2 = seg2->endpoint->x;
  242.  
  243.     if (s2>Cross_x || s1>Cross_x || Cross_x>e1 || Cross_x>e2) return false;
  244.  
  245.     //y-coordinate of the intersection
  246.     double Cross_y = (seg1->slope * Cross_x + seg1->yshift);
  247.     inter =  new SW_POINT(Cross_x,Cross_y);
  248.  
  249.     return true;
  250.   }
  251. }
  252.  
  253.  
  254. static pq_item Xinsert(seq_item i, SW_point p) 
  255. { 
  256.   return X_structure.insert(i,p);
  257. }
  258.  
  259. static SW_point Xdelete(pq_item i) 
  260. {
  261.   SW_point p = X_structure.inf(i);
  262.   X_structure.del_item(i);
  263.   return p;
  264. }
  265.  
  266.  
  267. static node New_Node(GRAPH<point,int>& G,double x, double y )
  268. { return G.new_node(point(x,y)); }
  269.  
  270. static void New_Edge(GRAPH<point,int>& G,node v, node w, SW_segment l )
  271. { if (get_orient(l)==0)
  272.        G.new_edge(v,w,get_color(l));
  273.   else G.new_edge(w,v,get_color(l));
  274. }
  275.  
  276.  
  277. void handle_vertical_segment(GRAPH<point,int>& SUB, SW_segment l)
  278. { 
  279.   SW_point p = new SW_POINT(get_x(get_startpoint(l)),get_y(get_startpoint(l)));
  280.   SW_point q = new SW_POINT(get_x(get_endpoint(l)),get_y(get_endpoint(l)));
  281.  
  282.   SW_point r = new SW_POINT(get_x(p)+double(1),get_y(p));
  283.   SW_point s = new SW_POINT(get_x(q)+double(1),get_y(q));
  284.  
  285.   SW_segment bot = new SW_SEGMENT(p,r,0,0);
  286.   SW_segment top = new SW_SEGMENT(q,s,0,0);
  287.  
  288.   seq_item bot_it = Y_structure.insert(bot,nil);
  289.   seq_item top_it = Y_structure.insert(top,nil);
  290.   seq_item sit;
  291.  
  292.   node u,v,w;
  293.   SW_segment seg;
  294.   
  295.  
  296.   for(sit=Y_structure.succ(bot_it); sit != top_it; sit=Y_structure.succ(sit))
  297.   { seg = Y_structure.key(sit);
  298.  
  299.     double cross_y = (get_slope(seg) * get_x(p) + get_yshift(seg));
  300.  
  301.     v = get_left_node(seg);
  302.     if (v==nil)
  303.     { w = New_Node(SUB,get_x(p),cross_y);
  304.       set_left_node(seg,w);
  305.      }
  306.     else
  307.     { double vx = SUB[v].xcoord();
  308.       if ( vx < get_x(p)-EPS2) 
  309.       { w = New_Node(SUB,get_x(p),cross_y);
  310.         New_Edge(SUB,v,w,seg);
  311.         set_left_node(seg,w);
  312.        }
  313.       else w = v;
  314.      }
  315.  
  316.     u = get_left_node(l);
  317.     if (u!=nil && u!=w) New_Edge(SUB,u,w,l);
  318.     set_left_node(l,w);
  319.  
  320.    }
  321.   
  322.   delete l;
  323.   delete top;
  324.   delete bot;
  325.     
  326.   Y_structure.del_item(bot_it);
  327.   Y_structure.del_item(top_it);
  328.  
  329.  }
  330.  
  331.  
  332. void Sweep_Segments(const list<segment>& L1, 
  333.                     const list<segment>& L2, GRAPH<point,int>& SUB)
  334. {
  335.   SW_point    p,inter;
  336.   SW_segment  seg, l,lsit,lpred,lsucc,lpredpred;
  337.   pq_item  pqit,pxmin;
  338.   seq_item sitmin,sit,sitpred,sitsucc,sitpredpred;
  339.  
  340.  
  341.   int count=1;
  342.  
  343.   //initialization of the X-structure
  344.  
  345.   segment s;
  346.  
  347.   forall(s,L1) 
  348.    { SW_point p = new SW_POINT(s.start());
  349.      SW_point q = new SW_POINT(s.end());
  350.      seg = new SW_SEGMENT(p,q,0,count++);
  351.      Xinsert(nil,get_startpoint(seg));
  352.    }
  353.  
  354.   count = -1;
  355.  
  356.   forall(s,L2) 
  357.    { SW_point p = new SW_POINT(s.start());
  358.      SW_point q = new SW_POINT(s.end());
  359.      seg = new SW_SEGMENT(p,q,1,count--);
  360.      Xinsert(nil,get_startpoint(seg));
  361.    }
  362.  
  363.  
  364.   count = 0;
  365.  
  366.   x_sweep = -MAXINT;
  367.   y_sweep = -MAXINT;
  368.  
  369.  
  370.   while(!X_structure.empty())
  371.   {
  372.     pxmin = X_structure.find_min();
  373.     p = X_structure.inf(pxmin);
  374.  
  375.     sitmin = X_structure.key(pxmin);
  376.  
  377.     Xdelete(pxmin);
  378.  
  379.  
  380.  
  381.     if (get_kind(p) == Leftend)
  382.  
  383.     //left endpoint
  384.     { 
  385.  
  386.       l = get_seg(p); 
  387.  
  388.       x_sweep = get_x(p);
  389.       y_sweep = get_y(p);
  390.  
  391.  
  392.       if (get_x(p) == get_x(get_endpoint(l)))
  393.         handle_vertical_segment(SUB,l);
  394.       else
  395.       {
  396.  
  397. /*
  398.       sit = Y_structure.lookup(l);
  399.       if (sit!=nil)  
  400.            error_handler(1,"plane sweep: sorry, overlapping segments");
  401. */
  402.  
  403.  
  404.       sit = Y_structure.insert(l,nil);
  405.  
  406.       Xinsert(sit,get_endpoint(l));
  407.  
  408.       sitpred = Y_structure.pred(sit);
  409.       sitsucc = Y_structure.succ(sit);
  410.  
  411.       if (sitpred != nil) 
  412.       { if ((pqit = Y_structure.inf(sitpred)) != nil)
  413.           delete Xdelete(pqit);
  414.  
  415.         lpred = Y_structure.key(sitpred);
  416.  
  417.         Y_structure.change_inf(sitpred,nil);
  418.  
  419.         if (intersection(lpred,l,inter))
  420.             Y_structure.change_inf(sitpred,Xinsert(sitpred,inter));
  421.       }
  422.  
  423.  
  424.       if (sitsucc != nil)
  425.       { lsucc = Y_structure.key(sitsucc);
  426.         if (intersection(lsucc,l,inter))
  427.            Y_structure.change_inf(sit,Xinsert(sit,inter));
  428.       }
  429.      } /* else if vertical */
  430.  
  431.     }
  432.     else if (get_kind(p) == Rightend)
  433.          //right endpoint
  434.          { 
  435.            x_sweep = get_x(p);
  436.            y_sweep = get_y(p);
  437.  
  438.            sit = sitmin;
  439.  
  440.            sitpred = Y_structure.pred(sit);
  441.            sitsucc = Y_structure.succ(sit);
  442.  
  443.            SW_segment seg = Y_structure.key(sit);
  444.  
  445.            Y_structure.del_item(sit);
  446.  
  447.            delete seg;
  448.  
  449.            if((sitpred != nil)&&(sitsucc != nil))
  450.            {
  451.              lpred = Y_structure.key(sitpred);
  452.              lsucc = Y_structure.key(sitsucc);
  453.              if (intersection(lsucc,lpred,inter))
  454.                 Y_structure.change_inf(sitpred,Xinsert(sitpred,inter));
  455.            }
  456.          }
  457.          else
  458.          /*point of intersection*/
  459.          { 
  460.            node w = New_Node(SUB,get_x(p),get_y(p));
  461.  
  462.            count++;
  463.  
  464.            /* Let L = list of all lines intersecting in p 
  465.  
  466.               we compute sit     = L.head();
  467.               and        sitpred = L.tail();
  468.  
  469.               by scanning the Y_structure in both directions 
  470.               starting at sitmin;
  471.  
  472.            */
  473.  
  474.            /* search for sitpred upwards from sitmin: */
  475.  
  476.            Y_structure.change_inf(sitmin,nil);
  477.  
  478.            sitpred = Y_structure.succ(sitmin);
  479.  
  480.  
  481.            while ((pqit=Y_structure.inf(sitpred)) != nil)
  482.            { SW_point q = X_structure.inf(pqit);
  483.              if (compare(p,q) != 0) break; 
  484.              X_structure.del_item(pqit);
  485.              Y_structure.change_inf(sitpred,nil);
  486.              sitpred = Y_structure.succ(sitpred);
  487.             }
  488.  
  489.  
  490.            /* search for sit downwards from sitmin: */
  491.  
  492.            sit = sitmin;
  493.  
  494.            seq_item sit1;
  495.            
  496.            while ((sit1=Y_structure.pred(sit)) != nil)
  497.            { pqit = Y_structure.inf(sit1);
  498.              if (pqit == nil) break;
  499.              SW_point q = X_structure.inf(pqit);
  500.              if (compare(p,q) != 0) break; 
  501.              X_structure.del_item(pqit);
  502.              Y_structure.change_inf(sit1,nil);
  503.              sit = sit1;
  504.             }
  505.  
  506.  
  507.  
  508.            // insert edges to p for all SW_segments in sit, ..., sitpred into SUB
  509.            // and set left node to w 
  510.  
  511.            lsit = Y_structure.key(sitpred);
  512.  
  513.            node v = get_left_node(lsit);
  514.            if (v!=nil) New_Edge(SUB,v,w,lsit);
  515.            set_left_node(lsit,w);
  516.  
  517.            for(sit1=sit; sit1!=sitpred; sit1 = Y_structure.succ(sit1))
  518.            { lsit = Y_structure.key(sit1);
  519.  
  520.              v = get_left_node(lsit);
  521.              if (v!=nil) New_Edge(SUB,v,w,lsit);
  522.              set_left_node(lsit,w);
  523.             }
  524.  
  525.            lsit = Y_structure.key(sit);
  526.            lpred=Y_structure.key(sitpred);
  527.            sitpredpred = Y_structure.pred(sit);
  528.            sitsucc=Y_structure.succ(sitpred);
  529.  
  530.  
  531.            if (sitpredpred != nil)
  532.             { 
  533.               lpredpred=Y_structure.key(sitpredpred);
  534.  
  535.               if ((pqit = Y_structure.inf(sitpredpred)) != nil)
  536.                 delete Xdelete(pqit);
  537.  
  538.               Y_structure.change_inf(sitpredpred,nil);
  539.  
  540.  
  541.               if (intersection(lpred,lpredpred,inter))
  542.                 Y_structure.change_inf(sitpredpred,
  543.                                        Xinsert(sitpredpred,inter));
  544.              }
  545.  
  546.  
  547.            if (sitsucc != nil)
  548.             {
  549.               lsucc=Y_structure.key(sitsucc);
  550.  
  551.               if ((pqit = Y_structure.inf(sitpred)) != nil)
  552.                 delete Xdelete(pqit);
  553.                  
  554.               Y_structure.change_inf(sitpred,nil);
  555.  
  556.               if (intersection(lsucc,lsit,inter))
  557.                   Y_structure.change_inf(sit,Xinsert(sit,inter));
  558.              }
  559.  
  560.  
  561. // reverse the subsequence sit, ... ,sitpred  in the Y_structure
  562.  
  563.            x_sweep = get_x(p);
  564.            y_sweep = get_y(p);
  565.  
  566.            Y_structure.reverse_items(sit,sitpred);
  567.  
  568.           delete p;
  569.  
  570.          } // intersection
  571.  
  572.   }
  573.  
  574.   X_structure.clear();
  575.  
  576. }
  577.  
  578. main()
  579. { 
  580.   int N = read_int("N = ");
  581.  
  582.   list<segment> seglist1,seglist2;
  583.   
  584.   while (N--) 
  585.   { double x1 = random(-1000,-100);
  586.     double y1 = random(-1000,1000);
  587.     double x2 = random(100,1000);
  588.     double y2 = random(-1000,1000);
  589.     seglist1.append(segment(x1,y1,x2,y2));
  590.    }
  591.  
  592.   GRAPH<point,int> SUB;
  593.  
  594.  
  595.   float T;
  596.  
  597.   T = used_time();
  598.   cout << "SWEEP_SEGMENTS (double)  :         ";
  599.   cout.flush(); 
  600.   SWEEP_SEGMENTS(seglist1,seglist2,SUB);
  601.   cout<< string(" # = %d time = %6.2f sec",SUB.number_of_nodes(), used_time(T));
  602.   newline;
  603.  
  604.   SUB.clear();
  605.  
  606.   T = used_time();
  607.   cout << "Sweep_Segments (leda-rational):    ";
  608.   cout.flush(); 
  609.   Sweep_Segments(seglist1,seglist2,SUB);
  610.   cout<< string(" # = %d time = %6.2f sec",SUB.number_of_nodes(), used_time(T));
  611.   newline;
  612.  
  613. return 0;
  614.  
  615. }
  616.